




#include"sl_auto_different.h"


// 这里只实现反向传播自动计算微分的框架，具体的计算求导代码还得自己编写。

// 这个计算图还真不能硬套这里的 stl_graph 的函数。这两者差别太大了。这里的节点得自己定义父子节点。


#include<math.h>




typedef struct ___variable_mat_
{
	char str[8];    // 字符串描述符。
	bool trainable; // 记录当前节点是否可训练，标签和样本不可训练，参数可训练，不可训练的节点就不用计算梯度了。
	bool is_calculateed;  // 标记前向计算结果的时候，该节点是否有被计算。

	mat operand[2]; // 0,1 分别是左操作数和右操作数
	mat result;     // 

	/// <summary>
	/// 首先，这个函数最终的计算结果是一个标量。因此，有几个变量就有几个梯度，梯度矩阵尺寸与操作数的矩阵一样。
	/// 对最终结果的梯度= pfinal/pgrad = pfinal/presult * presult/pgrad。 
	/// presult/pgrad一般是更高阶的张量，因此这里只保留pfinal/pgrad，presult/pgrad放到calculate里由用户自己实现。
	/// </summary>
	mat grad[2];   
	// 最终结果对当前节点的result的梯度= pfinal/presult,有这个梯度，就可以计算最终结果对操作数的梯度。
	mat grad_result;  
	// Tresult_or_Fjacobi == 1,前向计算result；Tresult_or_Fjacobi == 0，反向计算grad[2]
	void (*calculate)(bool Tresult_or_Fjacobi, struct ___variable_mat_* vf);
	// 两个操作数节点的指针，方便进行前向后向的计算。
	struct ___variable_mat_* operand_var[2];  

	// result,grad[2],grad_result,在创建节点的时候分配内存。
	
}varmat;

void varmat_allocate(varmat* vm,int g0_r,int g0_c, int g1_r, int g1_c, int r_r, int r_c, int gr_r, int gr_c)
{
	// result,grad[2],在创建节点的时候分配内存。
	mat_allocate(&vm->result, r_r, r_c);
	mat_allocate(&vm->grad[0], g0_r, g0_c);
	mat_allocate(&vm->grad[1], g1_r, g1_c);
//	mat_allocate(&vm->grad_result, gr_r, gr_c);
}

void varmat_deallocate(varmat* vm)
{
	mat_deallocate(&vm->result);
	mat_deallocate(&vm->grad[0]);
	mat_deallocate(&vm->grad[1]);
//	mat_deallocate(&vm->grad_result);
}


// 激活函数
void varopt_ReLU(bool Tresult_or_Fjacobi, varmat* vf)
{
	////ReLU函数   y = Max(a * x, x);  一般设 a = 0.01
#define ReLUalpha 0.01f  
	// 计算
	if (Tresult_or_Fjacobi) {

		real_t** dat = vf->operand[0].data;
		real_t** rdat = vf->result.data;
		int r = vf->operand[0].rows;
		int c = vf->operand[0].cols;
		for (int i = 0; i < r; i++) {
			for (int j = 0; j < c; j++) {
				rdat[i][j] = dat[i][j] > 0 ? dat[i][j] : ReLUalpha * dat[i][j];
			}
		}
	}
	else
	{
		real_t** dat = vf->operand[0].data;
		real_t** grad = vf->grad[0].data;
		real_t** rgrad = vf->grad_result.data;
		int r = vf->operand[0].rows;
		int c = vf->operand[0].cols;
		for (int i = 0; i < r; i++) {
			for (int j = 0; j < c; j++) {
				grad[i][j] += rgrad[i][j] * (dat[i][j] > 0 ? 1 : ReLUalpha);
			}
		}
	}
	vf->is_calculateed = true;
}
// 范数，这个节点一定是最后一个节点，输出operand[0]所有元素的平方和。
void varopt_norm(bool Tresult_or_Fjacobi, varmat* vf)
{
	// 计算 结果必然是个标量。
	if (Tresult_or_Fjacobi) {
		vf->result.data[0][0] = mat_norm(vf->operand[0]);
	}
	else
	{
		// 范数节点，一定是最后一个节点。
		real_t rst_grad = 1;  // vf->grad_result.data[0][0];

		real_t** dat = vf->operand[0].data;
		real_t** grad = vf->grad[0].data;
		int r = vf->operand[0].rows;
		int c = vf->operand[0].cols;
		rst_grad *= 2;
		for (int i = 0; i < r; i++) {
			for (int j = 0; j < c; j++) {
				grad[i][j] += rst_grad * dat[i][j];
			}
		}
	}
	vf->is_calculateed = true;
}
// 矩阵加法
void varopt_add(bool Tresult_or_Fjacobi, varmat* vf)
{
	// 计算
	if (Tresult_or_Fjacobi) {
		mat_add(vf->operand[0], vf->operand[1], vf->result);
	}
	else
	{
// 用指标求和的方式表述矩阵乘法偏导数：
// Aij + Bij = Cij  pC/pA = 1 ,  d/dAij = (d/Cij) *  pC/pA = (d/Cij) * 1 
//                  pC/pB = 1 ,  d/dBij = (d/Cij) *  pC/pB = (d/Cij) * 1  
		mat_add(vf->grad[0], vf->grad_result, vf->grad[0]);  // 累加梯度
		mat_add(vf->grad[1], vf->grad_result, vf->grad[1]);  // 累加梯度
	}
	vf->is_calculateed = true;
}
// 矩阵减法
void varopt_sub(bool Tresult_or_Fjacobi, varmat* vf)
{
	// 计算
	if (Tresult_or_Fjacobi) {
		mat_sub(vf->operand[0], vf->operand[1], vf->result);
	}
	else
	{
		mat_add(vf->grad[0], vf->grad_result, vf->grad[0]);  // 累加梯度
		mat_sub(vf->grad[1], vf->grad_result, vf->grad[1]);  // 累加梯度
	}
	vf->is_calculateed = true;
}
// 矩阵乘法
void varopt_mul(bool Tresult_or_Fjacobi, varmat* vf)
{
	// 计算
	if (Tresult_or_Fjacobi) {
		mat_mul(vf->operand[0], vf->operand[1], vf->result);
	}
	else
	{
		// 计算梯度
// 用指标求和的方式表述矩阵乘法偏导数：
// Aij * Bjk = Cik  pC/pA = Bjk ,  d/dAij = (d/Cik) *  pC/pA = (d/Cik) * Bjk 
//                  pC/pB = Aij ,  d/dBjk = (d/Cik) *  pC/pB = (d/Cik) * Aij  
// 这里需要记录的偏导数是  pC/pA = Bjk 或者  pC/pB = Aij  还是一个矩阵。只是这个矩阵在计算梯度的时候需要留意		
		real_t** gradA = vf->grad[0].data;
		real_t** gradB = vf->grad[1].data;
		real_t** gradC = vf->grad_result.data;
		real_t** A = vf->operand[0].data;
		real_t** B = vf->operand[1].data;

		int i, j, k, I, J, K;
		I = vf->operand[0].rows;
		J = vf->operand[0].cols;
		K = vf->operand[1].cols; 
		// 下面开始累加梯度
		for (int i = 0; i < I; i++) {
			for (int j = 0; j < J; j++) {
				for (int k = 0; k < K; k++) {
					gradA[i][j] += gradC[i][k] * B[j][k];
					gradB[j][k] += gradC[i][k] * A[i][j];
				}
			}
		}
	}
	vf->is_calculateed = true;
}
//矩阵除法必须是方阵
void varopt_inv(bool Tresult_or_Fjacobi, varmat* vf)
{
	// 计算 result = 1 / operand[0] ;    grad[1]参与辅助运算。
	if (Tresult_or_Fjacobi) {

		mat_copyto(vf->operand[0], vf->grad[1]);
		mat_inversion(vf->grad[1], vf->result, Epsilon);
	}
	else
	{		
// 用指标求和的方式表述矩阵乘法偏导数： 设 A＾-1 = B ,即： Aij * Bjk = I = Dik , Dik=(i==k)?1:0 
// 	dAij * Bjk + Aij * dBjk = 0 >>  因为是指标求和，所以没法分离。
// dB/dA = -(A^-2) >> dB = -(A^-2)*dA   >>  dBij = -(A^-2)ik * dAkj 
// 设 dAkj = Akj * 0.0001,   dBij = -(A^-2)ik * Akj * 0.0001 = -(A^-1)im * (A^-1)mk * Akj * 0.0001 = -(A^-1)ij * 0.0001  // 很合理
// 设 标量 e 对Bij的导数是 de/dBij , 则 e 对 A的全微分 de =  de/dB * dB = de/dBij * dBij = - de/dBij * (A^-2)ik * dAkj  
// 	所以：  de / dAkj = - de/dBij * (A^-2)ik  = 

		mat_mul(vf->result, vf->result, vf->grad[1]);  // grad[1] 是 没有用到的，但初始化节点的时候有分配内存。
		real_t** gradA = vf->grad[0].data;
		real_t** A_2 = vf->grad[1].data;  // A^-2 , 雅可比矩阵。
		real_t** gradB = vf->grad_result.data;

		int i, j, k, I, J, K;
		I = vf->operand[0].rows;
		J = vf->operand[0].rows;
		K = vf->operand[0].rows;
		// 下面开始累加梯度
		for (int i = 0; i < I; i++) {
			for (int j = 0; j < J; j++) {
				for (int k = 0; k < K; k++) {
					gradA[k][j] -= gradB[i][j] * A_2[i][k];
				}
			}
		}
	}
	vf->is_calculateed = true;
}
//当前数据节点
void varopt_ths(bool Tresult_or_Fjacobi, varmat* vf)
{
	// 计算
	if (Tresult_or_Fjacobi) {
		mat_copyfrom(vf->result, vf->operand[0]);
	}
	else
	{
		mat_add(vf->grad[0], vf->grad_result, vf->grad[0]);  // 累加梯度
	}
	vf->is_calculateed = true;
}



/// <summary>
/// 初始化运算节点。初始化节点，必须先初始化底层节点。即左右操作数的节点。
/// </summary>
/// <param name="current_var">当前的节点</param>
/// <param name="str">说明字符串（7个字符）</param>
/// <param name="operate">操作函数</param>
/// <param name="l_ope">左操作数指针</param>
/// <param name="r_ope">有操作数指针</param>
/// <param name="result_rows">当前节点结果行数</param>
/// <param name="result_cols">当前节点结果列数</param>
/// <param name="trainable">该节点是否训练</param>
void varmat_init(varmat* current_var, char* str, void (*operate)(bool, varmat*),
	varmat* l_ope, varmat* r_ope, int result_rows, int result_cols,bool trainable) {

	// 思路：初始化，可以知道当前节点的result尺寸，左右操作数指针可以知道操作数的尺寸，
	// 进而可以为result和梯度矩阵分配内存。

	current_var->calculate = operate;
	for (int i = 0; i < 7 && str[i] != 0; i++)	current_var->str[i] = str[i];
	current_var->str[7] = 0;
	current_var->trainable = trainable;
	int g0_r = 1, g0_c = 1, g1_r = 1, g1_c = 1, r_r = 1, r_c = 1, gr_r = 1, gr_c = 1;
	bool is_error = false;
	if (operate == varopt_ths) {
		current_var->operand_var[0] = NULL;
		current_var->operand_var[1] = NULL;
		g0_r = result_rows; g0_c = result_cols;
		g1_r = 1; g1_c = 1;
		gr_r = result_rows;  gr_c = result_cols;
		r_r = gr_r;  r_c = gr_c;
	}
	if (operate == varopt_add || operate == varopt_sub) {
		current_var->operand_var[0] = l_ope;
		current_var->operand_var[1] = r_ope;
		g0_r = l_ope->result.rows; g0_c = l_ope->result.cols;
		g1_r = r_ope->result.rows; g1_c = r_ope->result.cols;
		gr_r = result_rows;  gr_c = result_cols;
		r_r = gr_r;  r_c = gr_c;
		if (g0_r != g1_r || g0_c != g1_c || g1_r != r_r || g1_c != r_c)  is_error = true;
	}
	if (operate == varopt_mul) {
		current_var->operand_var[0] = l_ope;
		current_var->operand_var[1] = r_ope;
		g0_r = l_ope->result.rows; g0_c = l_ope->result.cols;
		g1_r = r_ope->result.rows; g1_c = r_ope->result.cols;
		gr_r = result_rows;  gr_c = result_cols;
		r_r = gr_r;  r_c = gr_c;
		if (g0_r != r_r || g0_c != g1_r || g1_c != r_c)   is_error = true;
	}
	if (operate == varopt_inv) {
		current_var->operand_var[0] = l_ope;
		current_var->operand_var[1] = NULL;
		g0_r = l_ope->result.rows; g0_c = l_ope->result.cols;
		g1_r = l_ope->result.rows; g1_c = l_ope->result.cols;
		gr_r = result_rows;  gr_c = result_cols;
		r_r = gr_r;  r_c = gr_c;
		if (g0_r != g0_c || r_r != r_c || r_r != g0_r)   is_error = true;
	}
	if (operate == varopt_ReLU) {
		current_var->operand_var[0] = l_ope;
		current_var->operand_var[1] = NULL;
		g0_r = l_ope->result.rows; g0_c = l_ope->result.cols;
		g1_r = 1; g1_c = 1;
		gr_r = result_rows;  gr_c = result_cols;
		r_r = gr_r;  r_c = gr_c;
		if (g0_r != r_r || r_c != g0_c)                 is_error = true;
	}
	if (operate == varopt_norm) {
		current_var->operand_var[0] = l_ope;
		current_var->operand_var[1] = NULL;
		g0_r = l_ope->result.rows; g0_c = l_ope->result.cols;
		g1_r = 1; g1_c = 1;
		gr_r = result_rows;  gr_c = result_cols;
		r_r = gr_r;  r_c = gr_c;
		if (result_rows != 1 || result_cols != 1)    is_error = true;
		
	}


	if (is_error) {
		printfs("\r\n  data size error !!  varmat: %s  init failed !  ", str);
		return;
	}


	// result,grad[2],grad_result,在创建节点的时候分配内存。
	varmat_allocate(current_var, g0_r, g0_c, g1_r, g1_c, r_r, r_c, gr_r, gr_c);
	l_ope = current_var->operand_var[0];
	r_ope = current_var->operand_var[1];
	if (l_ope != NULL) {
		current_var->operand[0] = l_ope->result;  // 操作数从子节点获取。 
		l_ope->grad_result = current_var->grad[0]; // 子节点的结果梯度从当前节点获取
	}
	if (r_ope != NULL) {
		current_var->operand[1] = r_ope->result;  // 操作数从子节点获取。 
		r_ope->grad_result = current_var->grad[1]; // 子节点的结果梯度从当前节点获取
	}

}

void varmat_deinit(varmat* current_var)
{
	varmat_deallocate(current_var);
}


void print_varmat(varmat* vf)
{
	printfs("\r\n %s  ", vf->str);
	printfs(" \r\n result:  ");
	mat_print(vf->result);

	printfs(" \r\n resule grad :  ");
	mat_print(vf->grad_result);
	printfs(" \r\n left operand:  ");
	mat_print(vf->operand[0]);
	printfs(" \r\n right operand:  ");
	mat_print(vf->operand[1]);
	printfs(" \r\n left grad:  ");
	mat_print(vf->grad[0]);
	printfs(" \r\n right grad:  ");
	mat_print(vf->grad[1]);

}

// 清空计算标志位
void varmat_forward_init(varmat* vf)
{
	// 如果操作数不空
	if (vf->operand_var[0] != NULL) {
		varmat_forward_init(vf->operand_var[0]);
	}
	if (vf->operand_var[1] != NULL) {
		varmat_forward_init(vf->operand_var[1]);
	}
	vf->is_calculateed = false;
}
// 计算节点结果。
void varmat_forward(varmat* vf)
{
	// 如果操作数不空，同时没有计算，先计算操作数
	if (vf->operand_var[0] != NULL) {
		if (vf->operand_var[0]->is_calculateed == false)
			varmat_forward(vf->operand_var[0]);
	//	vf->operand[0] = vf->operand_var[0]->result;
	}
	if (vf->operand_var[1] != NULL) {
		if (vf->operand_var[1]->is_calculateed == false)
			varmat_forward(vf->operand_var[1]);
	//	vf->operand[1] = vf->operand_var[1]->result;
	}
	// 计算当前节点。
	vf->calculate(true, vf);
//	print_varmat(vf);
}

// 清空计算标志位
void  varmat_backward_init(varmat* vf)
{
	// 如果操作数不空
	if (vf->operand_var[0] != NULL) {
		varmat_backward_init(vf->operand_var[0]);
	}
	if (vf->operand_var[1] != NULL) {
		varmat_backward_init(vf->operand_var[1]);
	}
	mat_init(vf->grad[0], 0);
	mat_init(vf->grad[1], 0);
}
// 计算反向计算节点梯度。
void varmat_backward(varmat* vf)
{
	vf->calculate(false, vf);
//	print_varmat(vf);

	// 如果操作数不空，回溯计算各个参数的梯度。
	if (vf->operand_var[0] != NULL && vf->operand_var[0]->trainable) {
	//	varmat_backward(vf->operand_var[0], vf->grad[0]);		
		varmat_backward(vf->operand_var[0]);		
	}

	if (vf->operand_var[1] != NULL && vf->operand_var[1]->trainable) {
	//	varmat_backward(vf->operand_var[1], vf->grad[1]);
		varmat_backward(vf->operand_var[1]);
	}

}



// 矩阵的自动微分也测试成功！！！！ 一旦有了思路，编码就成纯体力活了。
void test_auto_different()
{
	int size = 10;
	mat dataset, lable;   //  dataset 是 size*size的矩阵，每一个行向量表示一个数据，lable是单位矩阵。
	// 把dataset 看作是10个行矢量的数据集，有10个分类，用softmax分类函数，有10个输出。其中对应的lable行矢量就是输出的正确矢量	
	mat mA, mB, mC, mA2, mB2, mC2;
	
	mat_allocate(&dataset, 1, size);
	mat_allocate(&lable, 1, size);
	//mat_allocate(&dataset, size, size);
	//mat_allocate(&lable, size, size);
	mat_allocate(&mA, size, 4);
	mat_allocate(&mB, 4, size);
	mat_allocate(&mC, 1, 4);
	mat_allocate(&mA2, size, 4);
	mat_allocate(&mB2, 4, size);
	mat_allocate(&mC2, 1, 4);
	mat_init(lable, 0);  lable.data[0][0] = 1;
	mat_randomize(dataset, 23, -10, 10);
	mat_randomize(mA, 234, -1, 1);
	mat_randomize(mB, 345, -1, 1);
	mat_randomize(mC, 455, -1, 1);
	// 模型：min || ReLU(x*A + C) * B - lable ||  ,x是1*10，lable是1*10， A是10*4，B是4*10,C是1*4，
	varmat x, y, A, B, C, xA, xAC, RxAC, RxACB, RxACB_y, error;
	varmat_init(&x, "x", varopt_ths, NULL, NULL, 1, size, false); x.operand[0] = dataset;
	varmat_init(&y, "y", varopt_ths, NULL, NULL, 1, size, false); y.operand[0] = lable;
	varmat_init(&A, "A", varopt_ths, NULL, NULL, size, 4, true); A.operand[0] = mA;
	varmat_init(&B, "B", varopt_ths, NULL, NULL, 4, size, true); B.operand[0] = mB;
	varmat_init(&C, "C", varopt_ths, NULL, NULL, 1, 4, true); C.operand[0] = mC;
	varmat_init(&xA, "xA", varopt_mul, &x, &A, 1, 4, true);
	varmat_init(&xAC, "xAC", varopt_add, &xA, &C, 1, 4, true);
	varmat_init(&RxAC, "RxAC", varopt_ReLU, &xAC,NULL, 1, 4, true);
	varmat_init(&RxACB, "RxACB", varopt_mul, &RxAC, &B, 1, size, true);
	varmat_init(&RxACB_y, "RxACB_y", varopt_sub, &RxACB, &y, 1, size, true);
	varmat_init(&error, "error", varopt_norm, &RxACB_y, NULL, 1, 1, true);

	varmat_forward_init(&error);
	varmat_forward(&error);
	printfs(" \r\n  auto calculate: %f ", error.result.data[0][0]);

	varmat_backward_init(&error);
	varmat_backward(&error);
	printfs(" \r\n   auto calculate jacobi  ");
	mat_print(A.grad[0]);
	mat_print(B.grad[0]);
	mat_print(C.grad[0]);
	float record[9];
	record[0] = error.result.data[0][0];
	float delta = 0.000001f;// 0.0000001f;
	printfs(" \r\n   \r\n    grad A  ");
	for (int i = 1; i < 9; i++) {
		mat_mul_real(A.grad[0], delta, mA2);
		mat_sub(mA, mA2, mA);  // 这里mA一直在递减
		varmat_forward_init(&error);
		varmat_forward(&error);
		printfs(" \r\n %d   auto calculate: %f ", i, error.result.data[0][0]);
		record[i] = error.result.data[0][0];
	}
	for (int i = 1; i < 9; i++) {
		printfs(" \r\n delta ==  %f ", record[i] - record[i - 1]);
	}

	printfs(" \r\n   \r\n    grad B  ");
	for (int i = 1; i < 9; i++) {
		mat_mul_real(B.grad[0], delta, mB2);
		mat_sub(mB, mB2, mB);  // 这里一直在递减
		varmat_forward_init(&error);
		varmat_forward(&error);
		printfs(" \r\n %d   auto calculate: %f ", i, error.result.data[0][0]);
		record[i] = error.result.data[0][0];
	}
	for (int i = 1; i < 9; i++) {
		printfs(" \r\n delta ==  %f ", record[i] - record[i - 1]);
	}

	printfs(" \r\n   \r\n    grad C  ");
	for (int i = 1; i < 9; i++) {
		mat_mul_real(C.grad[0], delta, mC2);
		mat_sub(mC, mC2, mC);  // 这里mA一直在递减
		varmat_forward_init(&error);
		varmat_forward(&error);
		printfs(" \r\n %d   auto calculate: %f ", i, error.result.data[0][0]);
		record[i] = error.result.data[0][0];
	}
	for (int i = 1; i < 9; i++) {
		printfs(" \r\n delta ==  %f ", record[i] - record[i - 1]);
	}

	mat_deallocate(&dataset);
	mat_deallocate(&lable);
	mat_deallocate(&mA);
	mat_deallocate(&mB);
	mat_deallocate(&mC);
	mat_deallocate(&mA2);
	mat_deallocate(&mB2);
	mat_deallocate(&mC2);
	varmat_deinit(&x);
	varmat_deinit(&y);
	varmat_deinit(&A); 
	varmat_deinit(&B);
	varmat_deinit(&C); 
	varmat_deinit(&xA);
	varmat_deinit(&xAC);
	varmat_deinit(&RxAC);
	varmat_deinit(&RxACB);
	varmat_deinit(&RxACB_y);
	varmat_deinit(&error);

}


